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[i]  In  this  study,  we  present  an  aerosol  data  assimilation  system  destined  for  operational 
use  at  the  Fleet  Numerical  Meteorological  and  Oceanographic  Center  (FNMOC).  The 
system  is  an  aerosol  physics  version  of  the  Naval  Research  Laboratory  (NRL) 

Atmospheric  Variational  Data  Assimilation  System  (NAVDAS)  that  is  already  operational. 

The  purpose  of  this  new  system,  NAVDAS-Aerosol  Optical  Depth  (NAVDAS-AOD)  is  to 
improve  the  NRL  Aerosol  Analysis  and  Prediction  System  (NAAPS)’s  forecasting 
capability  by  assimilating  observational  data  sources  with  NAAPS  forecast  fields.  This 
will  allow  for  not  only  improved  aerosol  forecasting  but  also  for  dramatically  enhanced 
global  scale  research  capabilities  for  the  study  of  aerosol-meteorology  interaction. 

NAVDAS-AOD  assimilates  a  newly  developed  over-water  Moderate-Resolution  Imaging 
Spectroradiometers  (MODIS)  level  3  aerosol  product  with  NAAPS.  This  paper  is  the 
second  in  a  series  which  describes  NRL’s  program  to  realistically  monitor  global 
aerosol  distributions.  Here  we  explain  the  reasons  and  procedures  for  constructing  the 
over-water  level  3  MODIS  aerosol  product,  describe  the  theoretical  basis  for 
NAVDAS-AOD,  and  provide  a  thorough  statistical  error  analysis  for  both  the  MODIS 
observations  and  the  NAAPS  model  background  fields  that  are  critical  to  aerosol 
data  assimilation.  Using  5  months  of  analysis,  our  study  shows  that  by  carefully  screening 
over-water  satellite  observations  to  ensure  only  the  best  quality  data  are  used  in  the 
aerosol  assimilation  process,  the  NAVDAS-AOD  can  significantly  improve  the  NAAPS 
global  aerosol  optical  depth  analysis  as  well  as  improve  the  aerosol  forecast  skill. 

Citation:  Zhang,  J.,  J.  S.  Reid,  D.  L.  Westphal,  N.  L.  Baker,  and  E.  J.  Hyer  (2008),  A  system  for  operational  aerosol  optical  depth 
data  assimilation  over  global  oceans,  7.  Geophys.  Res.,  113,  D10208,  doi:10.l029/2007JD009065. 


1.  Introduction 

[2]  The  relatively  recent  use  of  multichannel,  multiangle 
observations  from  the  Moderate-Resolution  Imaging  Spec¬ 
troradiometers  (MODIS)  and  Multiangle  Imaging  Spectror- 
adiometer  (MISR)  instmments  has  dramatically  improved 
spatial  and  spectral  observations  of  aerosol  optical  depth 
(AOD).  Recent  developments  also  include  algorithmic 
advancements  (including  more  realistic  optical  models), 
better  aerosol  detection,  and  improvements  in  cloud  screen¬ 
ing  capabilities  [e.g.,  Kaufman  et  ai^  2002].  The  accuracy 
of  satellite  aerosol  optical  depth  retrievals  and  screening  has 
improved  to  a  point  where  numerical  aerosol  forecasting 
can  also  benefit  from  satellite  observations  through  the  data 
assimilation  process  [Zhang  and  Reid,  2006]. 

[3]  The  idea  of  assimilating  satellite  observations  into 
numerical  models  is  now  commonplace.  Some  successful 
applications  include  assimilation  of  a  satellite  wind  product 
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[Tomassini  et  al.,  1999],  satellite-retrieved  surface  temper¬ 
ature  [Dabberdt  and  Schlatter.,  1996],  and  direct  assimila¬ 
tion  of  satellite  radiances  [McNally  and  Vesperini.,  1996]. 
However,  only  a  few  attempts  have  been  made  to  assimilate 
satellite  aerosol  products  into  numerical  models  using  an 
optimal  interpolation  (01)  technique,  mainly  through  retro¬ 
spective  studies  or  limited  regional  studies  [e.g.,  Yu  et  ai, 
2004;  Collins  et  ai,  2001,  Rasch  et  aL,  2001]  and  most 
recently  through  radiance  assimilation  [Weaver  et  al..,  2007]. 
No  attempt  has  ever  been  made  to  assimilate  satellite 
aerosol  optical  depth  products  in  a  truly  operational  fashion 
globally.  This  is  largely  because  observations  and  numerical 
models  have  only  recently  improved  to  the  point  where 
serious  attempts  at  satellite  aerosol  assimilation  can  begin. 
Successful  data  assimilation  requires  that  the  observations 
and  the  forecast  model  are  sufficiently  accurate  and  have 
roughly  similar  quality.  Moreover,  operational  aerosol  data 
assimilation  requires  processed  data  to  be  accessible  within 
12  h  of  the  satellite  overpass  time. 

[4]  Recently,  Weaver  et  al.  [2007]  demonstrated  the 
concept  of  directly  assimilating  MODIS  radiances  in  into 
the  Goddard  Chemistry  and  Aerosol  Radiation  Transport 
model  (GOCART)  using  an  offline  assimilation  system. 
However,  given  the  operational  constraints  of  a  global 
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forecasting  system,  here  we  choose  to  use  the  conventional 
aerosol  products  that  arc  widely  accepted  through  the 
community.  We  optimistically  assume  that  the  algorithms 
are  accurate  and  that  instrument-specific  biases  have  been 
eliminated  by  the  developers  or  at  least  accounted  for  in  the 
reported  error  statistics  and  described  in  publications.  By 
using  aerosol  products  such  as  operational  MODIS  aerosol 
product,  we  are  taking  advantage  of  years  of  work  by  many 
different  groups  to  develop  and  validate  these  products  [e.g., 
Renter  et  al.^  2005;  Hsu  et  al.^  2006;  Husar  et  ai,  1997; 
Kahn  et  al.^  2005].  An  aerosol  assimilation  system  that  uses 
aerosol  optical  depth  (AOD)  instead  of  radiances  can  be 
easily  expanded  to  include  observations  from  other  existing 
or  future  sensors;  again,  provided  the  error  statistics  are 
known.  Validation  of  the  AOD  observations,  forecasts,  and 
analyses  is  more  easily  performed  since  AOD  is  a  com¬ 
monly  measured  quantity  (c.g..  Aerosol  Robotic  Network 
(AERONET)  global  network  of  Sun  photometers  [Holben  et 
al.,  1998]).  Last,  the  success  of  radiance  assimilation 
depends  on  the  accuracy  of  forecasted  temperature  and 
moisture  fields.  Assimilation  of  AOD  avoids  this  immediate 
dependency. 

[5]  In  any  data  assimilation  process,  high  data  quality  is 
required,  and  biases  and  noises  in  the  observational  data 
must  be  carefully  examined.  Large  uncertainties  exist  in 
aerosol  products,  especially  for  over  land  retrievals  where 
bright  heterogeneous  surface  characteristics  pose  a  chal¬ 
lenging  problem  [Levy  et  al.,  2007].  Studies  have  also 
shown  that  even  with  the  relatively  accurate  aerosol  retriev¬ 
als  over  oceans,  extra  quality  assurance  steps  need  to  be 
taken  before  implementation  of  the  aerosol  product  into 
models  [Zhang  and  Reid,  2006].  A  careful  data  screening 
process  is  needed  to  remove  noisy  data,  correct  biases,  and 
ensure  only  the  best  quality  data  are  used  in  an  aerosol  data 
assimilation  scheme. 

[6]  Currently,  there  are  three  aerosol  products  that  are 
created  in  a  near  real  time  mode  and  could  be  used  in  the 
data  assimilation  process:  MODIS,  AVHRR,  and  TOMS/ 
OMI.  Among  the  three  products,  we  chose  the  Aqua  and 
Terra  MODIS  aerosol  products  due  to  their  fine  spectral, 
spatial,  and  temporal  resolutions.  Comparing  with  previous 
sensors  such  as  AVHRR,  the  MODIS  products  have  dem¬ 
onstrated  superior  performance  in  both  cloud  and  aerosol 
detections  [Renter  et  al,  2005;  Zhang  and  Reid,  2006]. 

[7]  This  paper  is  the  second  in  a  series  which  outlines 
advances  in  the  Naval  Research  Laboratory’s  (NRL’s) 
global  aerosol  data  assimilation  and  modeling  program. 
The  first  provided  a  thorough  analysis  of  correction  factors 
for  MODIS  over-ocean  aerosol  optical  depth  products  for 
use  in  data  assimilation  [Zhang  and  Reid,  2006].  Here  we 
present  one  of  the  world’s  first  aerosol  data  assimilation 
systems  suitable  for  operation  use  that  assimilates  the 
quality  controlled  (QC)  and  quality  assured  (QA)  MODIS 
over-water  aerosol  product  described  in  that  previous  paper 
with  NAAPS  forecast  fields.  For  the  data  assimilation 
package,  we  adopted  and  modified  the  NRL  Atmospheric 
Variational  Data  Assimilation  System  (NAVDAS),  which 
has  been  used  operationally  for  assimilation  of  conventional 
and  satellite-based  observations  [Baker  et  ai,  2005].  In 
contrast  with  previous  aerosol  assimilation  efforts  [Collins 
et  ai,  2001;  Yu  et  ai,  2004],  we  use  a  2-D  variational 
approach.  In  comparison  to  the  01  technique,  the  2-D,  3-D, 


and  4-D  variational  techniques  arc  more  suitable  for  oper¬ 
ational  use  since  local  data  selection  is  unnecessary  and  a 
more  powerful  statistical  error  analysis  can  be  applied. 

[8]  In  this  manuscript,  we  present  an  outline  of  the 
NAVDAS  Aerosol  Optical  Depth  package  (NAVDAS- 
AOD).  Data  screening  processes  for  the  operational  MODIS 
over- water  aerosol  product  are  briefly  summarized.  Focus 
points  include  a  study  of  the  error  covariance  matrix  of  the 
NAAPS  model  (background)  field,  an  evaluation  of  the 
performance  of  5  months  of  MODIS  aerosol  assimilation 
using  NAVDAS-AOD,  and  the  impact  on  the  first  48  h  of 
NAAPS  aerosol  forecasts. 

2.  Satellite  Data 

[9]  In  NAVDAS-AOD,  Terra  (MOD04)  and  Aqua 
(MYD04)  MODIS  level  II  aerosol  products  are  collected 
from  the  NOAA  NESDIS  Near  Real  Time  Processing  Effort 
(NRTPE)  which  is  produced  with  a  latency  of  3  h  or  less  in 
most  cases.  NRTPE  is  a  joint  effort  from  NASA,  NOAA, 
Air  Force  Weather  Agency,  the  Naval  Research  Laboratory, 
and  the  Naval  Oceanographic  Office.  With  its  reliable  data 
link  and  short  data  latency,  NRTPE  has  proven  to  be  one  of 
the  best  data  sources  for  operational  near-real  time  data.  For 
the  study  presented  here,  NAVDAS-AOD  is  tested  with 
5  months  of  data  collected  from  January  to  May  2006. 

[10]  MODIS  has  a  total  of  36  spectral  channels  ranging 
from  near-ultraviolet  to  infrared.  Over  global  oceans,  an 
optimal  minimization  technique  is  applied  to  six  spectral 
channels  (ranging  from  0.55  to  1.2  ^m)  for  retrieving  aerosol 
optical  depth  (r)  and  aerosol  size  parameter  [Renter  et  ai, 
2005].  Although  MODIS  can  provide  aerosol  retrievals  at 
multiple  channels,  the  aerosol  r  in  this  paper  refers  to  r  at 
0.55  fim.  Although  only  one  channel  information  is  used  in 
this  study,  the  utilization  of  the  spectral  signatures  from 
multiple  channels  is  important  and  is  subject  to  future  studies. 
Renter  et  ai  [2005]  gave  an  overall  accuracy  for  over- water 
aerosol  optical  depth  of  0.03  ±  0.05  r.  However,  Zhang  and 
Reid  [2006]  showed  that  systematic  biases  exist  in  MODIS 
over-water  aerosol  retrievals.  These  are  related  to  cloud 
fraction  and  contamination,  aerosol  type,  and  near-surface 
wind  speed  impacting  the  lower  boundary  condition. 
Empirical  corrections  were  proposed  and  data  screening 
strategies  were  applied  to  remove  outliers.  In  this  study,  we 
adopted  the  empirical  corrections  and  data  screening  pro¬ 
cesses  described  by  Zhang  and  Reid  [2006]  where  a  total  of 
40%  data  were  lost  but  70%  of  the  significant  outliers  were 
removed.  This  strict  QA  procedure  is  necessary  to  remove 
these  outliers  that  would  adversely  impact  the  model  and 
propagate  through  the  system,  particularly  in  the  middle-  to 
high-latitude  oceans. 

[11]  After  the  empirical  corrections  and  data  screening 
procedures,  MODIS  r  data  is  binned  into  1®  x  1°  latitude/ 
longitude  grids  to  form  a  6  hourly,  data  assimilation  quality. 
Level  3  product  (e.g.,  at  NAAPS  time  of  t  =  6  h,  we  used 
both  Terra  and  Aqua  MODIS  data  from  t  =  3  h  to  t  =  9  h).  In 
this  step,  we  further  exclude  data  entries  with  MODIS  cloud 
fraction  larger  than  30%.  We  also  add  an  additional  QC 
procedure  by  removing  bins  that  are  isolated  from  other 
continuous  aerosol  features  and  bins  with  less  than  five 
valid  data  entries.  The  additional  QC  procedures  are  includ¬ 
ed  to  minimize  erroneous  aerosol  features  in  MODIS 
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aerosol  product  that  are  introduced  by  cloud  artifacts 
[Zhang  et  al.,  2005;  Zhang  and  Reid,  2006]. 

[12]  The  QA  and  QC  procedures  described  in  this  paper, 
however,  are  different  from  traditional  data  assimilation 
approaches  where  noises  in  observations  are  removed  based 
on  innovation  (difference  in  observations  and  forecasts) 
checks  [e.g..  Dee  et  ai,  2001].  We  argue  that  external  QA 
and  QC  checks  are  necessary  for  aerosol  data  assimilation 
because  of  the  following  reasons: 

[13]  1 .  Innovation  checks  could  not  flag  cases  when  both 
observations  (O)  and  forecasts  (F)  have  large  biases,  yet  O 
minus  F  values  are  small. 

[14]  2.  The  efficiency  of  observations  noise  removal 
could  be  sensitive  to  arbitrary  thresholds  used  in  the 
innovation  checks. 

[15]  3.  Large  innovation  values  could  be  valid  in  cases 
newly  formed  aerosol  plumes  with  significant  AOD  values 
which  models  fail  to  predict.  Examples  of  such  cases  are 
biomass  burning  events  with  fire  hot  spots  under  cloud 
covers.  Aerosol  forecasting  models  such  as  NAAPS  rely  on 
satellite  detected  hot  spots  to  define  their  emission  sources 
for  biomass  burning  aerosols.  Missing  satellite-detected  hot 
spots  will  cause  models  to  miss  a  possible  significant 
biomass  burning  episode.  However,  with  bias  correction, 
it  is  possible  to  include  such  cases,  as  will  be  discussed  in 
later  part  of  the  text  [Dee  and  da  Silva,  1999]. 

[16]  4.  By  performing  QA  with  this  method,  we  can 
generate  a  data  assimilation  quality  level  3  product  that  is 
applicable  to  any  model,  not  just  NAAPS/NAVDAS-AOD. 

3.  Description  of  the  NAAPS  Model 

[17]  With  its  transition  to  the  Fleet  Numerical  Meteorology 
and  Oceanography  Center  (FNMOC)  in  October  of  2006, 
NAAPS  is  the  U.S.  Navy’s  (and  the  world’s  first)  truly 
operational  global  aerosol  forecast  model.  NAAPS  produces 
6-d  forecasts  of  SO2,  sulfate,  dust,  biomass  burning  smoke 
and  sea  salt  mass  concentration  with  1  ®  x  1  °  resolution  on  30 
vertical  levels  [Witek  et  ai,  2007]. 

[is]  The  NAAPS  global  model  is  a  modified  form  of  a 
hemispheric  model  of  sulfate  aerosols  developed  by 
Christensen  [1997].  NAAPS  is  an  offline  model  that  utilizes 
the  meteorological  analysis  and  forecast  fields  from  the  Navy 
Operational  Global  Analysis  and  Prediction  System 
(NOGAPS)  [Hogan  and  Rosmond,  1991;  Hogan  and  Brody, 
1993].  Dynamical  fields  from  NOGAPS  are  remapped  to  1° 
resolution  from  the  0.5°  native  NOGAPS  resolution. 
Twice  daily,  the  NOGAPS  weather  forecast  model  provides 
dynamical  fields  to  NAAPS  at  6-h  intervals  for  the  6-d  forecast 
period.  Transport  is  calculated  using  a  3-D  semi-Lagrangian 
scheme  [Staniforth  and  Cote,  1991]  with  departure  points 
calculated  using  the  method  of  Ritchie  [1987]. 

[19]  Modifications  have  been  made  to  the  interpolation  of 
wind  and  concentration  fields  across  the  poles,  and  the 
interpolation  method  was  changed  from  a  third -order 
Lagrangian  to  fifth-order  Lagrangian.  Horizontal  and  ver¬ 
tical  diffusion  are  calculated  with  a  finite  element  scheme. 
The  vertical  diffusion  coefficient  parameterization  Kz  is 
based  on  the  Mon  in-Obukhov  similarity  theory  for  the 
surface  layer  using  the  input  NOGAPS  data.  The  Kz  profile 
is  extended  to  the  whole  boundary  layer  by  using  a  simple 
extrapolation  height  of  the  mixing  layer.  The  horizontal 


diffusion  coefficient  is  set  to  a  nominal  value  of  6  x  1 0"^  m^  s  “  * . 
The  condensation/precipitation  scheme  is  derived  from 
the  NOGAPS  atmospheric  model  and  is  further  described 
in  the  work  of  Hogan  and  Rosmond  [1991].  These  cloud 
profiles  are  used  to  calculate  the  wet  removal  and  reaction 
rates  of  the  various  species.  The  dry  deposition  velocity  is 
based  on  the  resistance  method  [Voldner  et  al.,  1986; 
Walcek  et  al.,  1986;  Slinn  and  Slinn,  1980],  where  the 
deposition  velocity  depends  on  the  turbulence  of  surface 
layer  and  surface  type  (ocean,  grassland,  etc.). 

[20]  Aerosol  source  functions  come  from  a  variety  of 
static  and  dynamical  algorithms.  Sulfur  dioxide  emission  is 
based  on  the  GEIA  inventory,  version  1  A,  for  the  year  1985 
with  a  seasonal  variation  and  two-level  vertical  distribution 
[Benkovitz  et  al.,  1996].  Natural  emissions  of  DMS  are 
immediately  converted  to  95%  sulfur  dioxide  and  5% 
sulfate.  The  gas-phase  chemistry  is  described  by  a  simple 
linear  reaction  rate,  which  depends  on  the  time  of  year  and 
latitude  [Christensen,  1997].  Dust  emission  occurs  when¬ 
ever  the  NOGAPS  friction  velocity  exceeds  a  threshold 
value  and  the  surface  moisture  is  less  than  0.3.  The 
threshold  friction  velocity  is  set  to  infinity  except  in  known 
dust-emission  areas  where  it  is  60  cm  s~^  [Westphal  et  al., 
1988].  These  areas  are  currently  defined  as  areas  covered  by 
the  erodible  land-use  types  used  in  the  USGS  Land  Cover 
Characteristics  Database  with  modification  guided  by  an 
analysis  of  TOMS  aerosol  index.  Smoke  emissions  are 
based  on  the  Fire  Locating  and  Modeling  of  Burning 
Emissions  (FLAMBE)  source  functions  based  on  near-real 
time  geostationary  and  MODIS  fire  products  [Reid  et  al., 
2004]. 

4,  Overview  of  the  Aerosol  Data  Assimilation 
Package 

[21]  NAVDAS-AOD  is  founded  on  the  3-D  variational 
analysis  option  of  the  NRL  Atmospheric  Variational  Data 
Assimilation  System  (NAVDAS)  [Daley  and  Barker,  2001]. 
NAVDAS  is  an  operational  three-dimensional  variational 
data  assimilation  suite  for  generating  atmospheric  state 
estimates  to  satisfy  a  variety  of  U.S.  Navy  needs.  NAVDAS 
is  formulated  in  observation  space.  The  preconditioned 
conjugate  gradient  method  [Daley  and  Barker,  2001]  is 
used  to  minimize  the  cost  function,  which  can  be  under¬ 
stood  as  the  process  to  minimize  the  analysis  error  variance 
in  linear  cases.  The  number  of  iterations  required  to  reach 
convergence  is  minimized  through  the  use  of  dual  block 
diagonal  preconditioners  with  Choleski  decomposition. 
Forward  operators  are  formulated  and  used  for  the  direct 
assimilation  of  Advanced  Television  Infrared  Observation 
Satellite  Operational  Vertical  Sounder  (ATOVS)  radiances 
and  SSM/I  wind  speeds  and  total  precipitable  water. 
NAVDAS  also  contains  a  comprehensive  diagnostic  suite, 
including  complete  observation  traceability,  Web-based 
observation  monitoring,  chi-square  monitoring  of  innova¬ 
tions,  the  adjoint  of  the  assimilation  system,  and  analysis 
error  variance  estimation. 

[22]  The  traditional  objective  in  data  assimilation  is  to 
create  a  new  analysis  field  by  correcting  the  background 
field  (i.e.,  the  current  model  stage)  using  observations,  and 
then  use  the  new  analysis  field  as  an  initial  condition  for 
future  model  runs.  In  the  case  of  the  assimilation  of  a  two- 
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dimensional  variable  such  as  optical  depth  (r),  we  degraded 
NAVDAS  from  3-D  to  2-D  variational  approach  with  three 
steps: 

[23]  1 .  Preprocessing:  a  process  to  convert  NAAPS  mass 
concentration  to  aerosol  optical  depth: 

TbA  =  )  +  £"bA  (  I  ) 

where,  TbA  is  the  background  (prior  forecast)  aerosol  optical 
depth  vector,  Cm  is  the  NAAPS  mass  concentration,  and 
H,n_T  is  the  forward  operator  that  represents  the  transforma¬ 
tions  from  NAAPS  mass  concentration  to  optical  depth  (see 
section  4.1  for  detail  discussions  of  H^  t)-  ^bA  is  the  error  in 
TbA  introduced  by  the  //^  r  operator. 

[24]  2.  Two-D  variational  assimilation  of  the  optical  depth 
field 

Tax  =  T,x  +  PiH^'  ^  /?] -'[Tax  -  H(t*;,)],  (2) 

where  TaA  is  the  analysis  optical  depth  vectors,  Tqa  is  the 
observation  optical  depth  vector,  and  H  is  the  observation 
operator  that  represents  any  necessary  spatial  and  temporal 
interpolations  from  background  to  the  observational  space. 
Pb  and  R  are  the  background  error  covariance  and  the 
observation  error  covariance  matrices,  respectively. 

[25]  The  analysis  field  (TaA)  can  be  regarded  as 
bacl^round  (TbA)  pius  a  correction  term  (ta,  also 
PftH^HP^H^+/?]“*[ToA-//('rbA)])»  where  the  correction 
term  is  the  difference  between  the  observation  and  back¬ 
ground  vectors  weighted  by  the  ratio  of  background  error 
covariance  matrix  to  total  error  covariance  matrix  in  the 
observational  space  (with  H  operator). 

[26]  3.  Postprocess:  a  process  to  convert  TaA  to  NAAPS 
mass  concentration: 

~  “I"  (3) 

where  //r  m  is  the  backAvard  operator  that  represents  the 
transformations  from  optical  depth  to  NAAPS  mass 


concentration,  is  the  error  in  Cm  introduced  by  the 
Hr  m  operator.  Details  of  the  //r_m  operator  will  be  discussed 
in  section  4.3.  Both  ^m  and  ^bA  could  be  transformed  as  part 
of  the  bias  term  of  TbA,  which  is  assumed  to  be  zero  for  this 
study.  The  background  bias  correction  will  be  discussed 
further  in  the  later  part  of  this  paper. 

[27]  The  NAAPS  prognostic  variable  is  the  3-D  aerosol 
mass  concentration.  In  this  study,  we  reduced  the  problem  to 
a  2-D  variational  assimilation  approach  simply  because 
the  operational  MODIS  aerosol  product  only  provides 
column  integrated  aerosol  optical  properties  (2-D).  With 
the  use  of  3-D  observations  of  aerosol  optical  properties 
from  CALIPSO  or  through  the  height  information  available 
from  some  gas  retrievals,  it  is  feasible  in  some  circum¬ 
stances  to  extend  equation  (2)  to  a  3-D  variational  approach 
in  future  studies. 

[28]  A  flowchart  of  the  NAVDAS-AOD  process  is  shown 
in  Figure  1.  Four  major  steps  are  included:  (1)  Convert 
NAAPS  mass  concentration  to  r  bA  (3-D  to  2-D  conver¬ 
sion);  (2)  run  NAVDAS  2-D  var  to  create  a  new  analysis 
field  TaA  frorn  r  bA  and  Tqa;  (3)  improve  the  NAAPS  mass 
concentration  field  using  TaA  (2-D  to  3-D  conversion);  and 
(4)  use  the  new  mass  concentration  field  as  an  initial 
condition  for  the  next  6-h  NAAPS  run. 

4.1.  NAAPS  Aerosol  Mass  to  Optical  Depth  Transfer 
Function  (//r_m  operator) 

[29]  In  the  first  step,  Tb  values  are  estimated  from  NAAPS 
3-D  mass  concentration  (g/m^)  using  a  bulk  model: 

TiA  =  5^o«(A)'C;/„(/J//,A)  (4) 

i 

where  for  each  species  /  and  wavelength  A,  TbA  is  the 
wavelength  dependent  background  aerosol  optical  depth, 
ae(A)  is  the  wavelength  dependent  mass  extinction 
efficiency  (m^  g”*),  and  Cm  is  the  aerosol  particle  mass 
concentration  (g  m”^)  and  yi(RH,  A)  is  the  aerosol 
hydroscopic  growth  factor  for  extinction.  For  this  first 
generation  model,  mass  extinction  cross  sections  (as  a 
function  of  wavelength)  for  dry  dust,  sulfate,  and  sea  salt  are 
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Figure  1.  Flowchart  of  the  NAVDAS-AOD. 


4  of  13 


D10208 


ZHANG  ET  AL.:  AEROSOL  DATA  ASSIMILATION 


D10208 


obtained  from  Optical  Properties  of  Aerosols  and  Clouds 
(OPAC)  [Hess  et  al^  1998].  The  parameters  used  in  this 
section  of  study  may  be  updated  when  new  observational 
parameters  are  available.  Smoke  properties  are  derived  from 
Reid  et  al  [2005].  The  aerosol  hygroscopic  growth  factor 
for  aerosol  absorption  coefficient  and  dust  scattering 
coefficient  are  assumed  to  be  zero.  The  aerosol  hygroscopic 
growth  factor  for  sulfate  and  smoke  scattering  coefficients 
are  estimated  based  on  Hanel  [1976]  and  Reid  et  al.  [2005], 
respectively.  The  aerosol  hygroscopic  growth  factor  for  sea 
salt  extinction  coefficient  is  estimated  from  Hegg  et  al. 
[2002]  and  Ming  and  Russell  [2001]. 

[30]  It  is  well  understood  that  that  aerosol  physical  and 
optical  properties  vary  regionally  and  are  highly  event 
based.  Therefore,  using  generalized  aerosol  optical  proper¬ 
ties  could  introduce  uncertainties  {e\,x)  in  the  estimated  Tb 
values.  But  for  the  purpose  of  radiation  studies,  aerosol 
optical  depth  and  extinction  are  the  variables  of  interest.  The 
use  of  a  bulk  model  in  this  way,  while  at  time  unphysical 
with  respect  to  aerosol  mass,  ensures  consistency  in  the 
atmospheric  radiation  fields.  Indeed,  the  area  of  mass 
extinction  efficiencies  and  extinction  to  mass  transfer  func¬ 
tion  is  an  area  of  intensive  research  in  the  community. 

4.2.  Creation  of  the  Optical  Depth  Analysis  Field 

[31]  In  the  second  step,  both  observation  (Tqa,  2-D)  and 
the  background  (xbA,  2-D)  vectors  developed  in  step  1  are 
interpolated  to  the  observational  domain,  and  a  new  analysis 
vector  (TbA  +  ta,  2-D)  is  created  by  correcting  background 
vector  from  observations  as  showed  in  equation  (2). 
To  estimate  ta,  both  background  and  observational  error 
covariance  matrices  are  required.  In  this  study,  background 
and  observational  error  covariance  are  defined  as  spatial 
error  covariance  and  are  assumed  to  be  independent 
from  meteorological  parameters  such  as  wind  speed.  The 
nature  of  the  observational  and  background  error  variance 
and  error  covariance  fields  will  be  discussed  in  detail  in 
section  5. 

4.3.  Reciprocal  Analysis  Optical  Depth  to  Mass 
Concentration  Transfer  Function  (H.r_m  Operator) 

[32]  The  third  step  in  the  optical  depth  data  assimilation 
process  is  to  apply  ta  to  update  the  NAAPS  mass  concen¬ 
tration  field.  While  it  is  relatively  straightforward  to  go  from 
3-D  to  2-D  as  in  section  4. 1 ,  without  additional  information 
the  converse  is  ill-posed.  To  correct  the  background  to  the 
analysis  field,  where  do  we  add  or  subtract  aerosol  particle 
mass  (in  vertical  and  compositional  dimensions)?  Further, 
how  do  we  know  the  differences  between  errors  in  mass  and 
hygroscopicity  (either  through  the  parameter  or  the 
NOGAPS  RH  field  itself)? 

[33]  Our  current  solution  is  to  keep  the  transfer  function 
straightforward  and  traceable.  If  f  a  is  less  than  zero,  which 
indicates  TbA  is  higher  than  the  ratio  of  {r^x  +  ^x)Hbx  is 
used  to  reduce  the  amplitude  of  mass  concentration  vertical 
profile  of  a  given  air  column;  i.e.,  the  concentration  in  each 
vertical  level  is  reduced  by  the  same  fractional  amount,  and 
same  for  each  aerosol  composition.  If  ta  =  0,  r^x  equals  to 
TbA,  ^nd  no  change  is  made  to  the  air  column.  If  rx  >  0,  TbA 
is  lower  than  TaA,  then  the  amplitude  of  the  vertical  profile 
of  aerosol  mass  is  uniformly  increased.  Given  NAAPS’  goal 
of  forecasting  large  visibility-reducing  events,  under  most 


circumstances  this  simple  method  is  sufficient.  Again,  with 
only  2-D  satellite  AOD  observations  available,  this  simple 
scaling  technique  is  a  valid  first  guess  to  redistribute  column 
integrated  AOD  vertically. 

[34]  Difficulties  arise  if  the  NAAPS  aerosol  source  func¬ 
tion  or  transport  meteorology  totally  misses  an  event  that  is 
observed  in  MODIS.  Newly  introduced  aerosol  features 
could  be  totally  different  from  the  model  predictions  in 
terms  of  aerosol  species  and  vertical  distributions.  One  of 
the  solutions  for  this  problem  is  to  adopt  a  background  bias 
correction  technique  [e.g..  Dee  and  da  Silva^  1998],  which 
is  the  subject  of  ongoing  research.  As  the  first  step  of  this 
study,  two  procedures  are  used  to  compensate  for  the 
problem:  (1)  for  TbA  >  0.10,  the  ratio  of  (TbA  +  ^xYr^x  is 
still  used  to  scale  up  the  mass  concentration  of  a  given  air 
column;  (2)  for  TbA  <  0.10,  a  seasonal  3-D  aerosol  clima¬ 
tology  created  using  3  years  (2004  2006)  of  NAAPS  data 
is  used.  The  NAAPS  3-D  aerosol  climatology  was  con¬ 
structed  for  dust,  smoke,  and  sulfate  aerosols,  for  broad 
regions  (60°  longitude  by  30°  latitude). 

4.4.  Forecast 

[35]  As  shown  in  Figure  1 ,  we  separated  NAAPS  model 
runs  into  two  modes:  assimilation  run  mode  and  forecast  run 
mode.  The  assimilation  run  mode  is  defined  as  the  NAAPS 
run  initialized  with  NAVDAS-AOD  analyses  every  6  h  in¬ 
terval.  The  free-running  forecast  or  nonassimilation  run 
mode  refers  to  the  time  step  when  NAAPS  model  runs  in 
a  forecast  mode  without  invoking  the  NAVDAS-AOD 
procedures  described  in  sections  4. 1-4.3.  Since  forecast 
accuracy  is  more  important  for  operational  models,  in  this 
study,  we  evaluated  the  impact  of  aerosol  data  assimilation 
not  only  to  NAAPS  assimilation  runs,  but  also  to  NAAPS 
forecasts. 

5.  Derivation  of  the  Error  and  Error  Covariance 
Matrices 

[36]  As  discussed  in  section  4,  both  background  (model 
space)  and  observational  error  covariance  matrices  are  crit¬ 
ical  to  the  aerosol  data  assimilation  process  and  arc  tightly 
connected  through  the  physics  of  measurement.  The  error 
covariance  matrices  are  composed  of  diagonal  terms  and 
nondiagonal  terms.  The  diagonal  terms  represent  error  var¬ 
iances,  and  the  nondiagonal  terms  represent  spatial  error 
covariances.  In  the  following  section,  we  show  our  estimates 
of  both  diagonal  and  nondiagonal  terms  of  the  observational 
error  covariance  matrix  in  section  5.1  and  show  the  NAAPS 
model  error  covariance  matrix  in  section  5.2. 

5.1.  Observational  Coverage  and  Errors 

[37]  In  this  study,  we  assume  that  observational  errors  are 
uncorrelated.  That  is,  the  error  in  one  1°  x  1°  MODIS 
level  3  aggregate  is  uncorrelated  to  the  ones  adjacent  to  it. 
In  reality,  this  is  certainly  not  the  case.  Because  the  MODIS 
algorithm  is  sensitive  to  aerosol  microphysics  (e.g.,  smoke 
versus  dust)  and  the  lower  boundary  condition  (glint,  white- 
capping),  any  large  aerosol  feature  will  have  correlated 
errors.  However,  in  the  QA  and  correction  procedures  of 
Zhang  and  Reid  [2006],  most  of  the  systematic  uncertainties 
have  been  corrected.  While  this  data  treatment  is  by  no 
means  perfect,  it  does  reduce  correlated  errors. 


5  of  13 


D10208 


ZHANG  ET  AL.:  AEROSOL  DATA  ASSIMILATION 


D10208 


[38]  The  reason  why  this  assumption  is  important  is  that  it 
will  zero  out  all  the  nondiagonal  elements  of  the  observa¬ 
tional  error  covariance  matrix,  which  leaves  only  diagonal 
elements  as  the  nonzero  terms.  The  diagonal  elements  are 
simply  observational  error  variances  and  are  estimated  using 
both  instrumental  error  variance  (erf)  and  sample  error 
variance  (<jf)  or  so-called  spatial  representative  error  vari¬ 
ance.  To  estimate  erf,  we  used  the  RMS  instrument  error 
variance  as  suggested  by  Zhang  and  Reid  [2006],  which  is 
similar  after  corrections  to  that  proposed  by  Renter  et  al. 
[2005]  for  T  <  0.6.  The  erf  represents  spatial  data  variation 
and  is  estimated  by  the  spatial  sample  variance  from  the 
averaging  of  MODIS  r  at  1°  x  1®  area.  The  total  observa¬ 
tional  error  variance  (o^)  is  simply  defined  as  equation  (5): 

(5) 

[39]  In  contrast  with  some  of  the  studies  that  estimate 
background  and  observational  error  variances  and/or  disen¬ 
tangle  the  individual  components  of  the  observation  errors 
using  innovation  statistics  [e.g..  Dee  and  da  Silva,  1998, 
1999],  we  used  AERONET  optical  depth  observations  that 
are  independent  from  both  model  and  satellite  observations. 
The  drawback  of  this  approach  is  that  AERONET  observa¬ 
tions  arc  point  source  observations  and  may  not  represent 
the  bin-averaged  AOD  values  in  some  cases.  However, 
AERONET  data,  which  is  not  included  in  the  assimilation 
process,  provide  independent  evaluations  of  both  back¬ 
ground  and  satellite  error  variances. 

5.2.  Background  Error  and  Error  Covariance  Matrix 

[40]  The  model  background  error  covariance  matrix  is 
critical  to  the  aerosol  data  assimilation  as  it  determines  the 
impact  ranges  of  observations.  For  any  given  two  grid 
locations,  m  and  n,  where  for  the  1®  NAAPS  model  m 
and  n  ranging  from  1  to  360  x  1 80,  the  background  error 
covariance  between  the  two  grid  points  (Pt^)  can  be 
estimated  as: 

P","  (6) 

where  and  Si  are  error  variances  at  location  m  and  n, 
respectively  (also  the  diagonal  terms  of  the  error  covariance 
matrix),  and  Q  is  the  error  correlation  between  point  m  and 
n.  The  NAAPS  error  variance  (5^)  is  estimated  as  a  function 
of  NAAPS  r,  based  on  half  year  (January- June  2006)  of 
analysis  using  NAAPS  and  AERONET  data.  To  estimate  S^, 
we  first  computed  model  error  variances  for  every  0.1 
NAAPS  T  interval  (for  r  <  1.0)  and  then  applied  a  linear 
regression  through  these  error  variance  values.  The  overall 
square  root  of  NAAPS  model  error  variance  is  estimated 
to  be  0.20  +  QAt  for  NAAPS  nonassimilation  runs  and  is 
0. 1 5  +  0.3r  for  NAAPS  6  h  forecast  with  assimilation. 

[41]  Give  the  fact  that  Q  contains  64,800  x  64,800  array 
elements  and  is  computationally  expensive  to  estimate  each 
single  term,  we  estimated  Q  using  the  second  order  autor¬ 
egressive  (SOAR)  function  [Daley  and  Barker,  2001],  as  the 
horizontal  correlation  model  and  is  shown  in  equation  (7): 

Cb(m,n)  =  (1  -r-  Rmn/ L)exp(-R„^„/L)  (7) 


where  is  the  great  circle  distance  and  L  is  the  global 
averaged  horizontal  error  correlation  length.  A  value  for  L 
of  200  km  is  used  in  this  study.  We  validated  our  choice  of  L 
in  section  6.3  by  studying  the  spatial  correlation  of 
observation  minus  forecast. 

6.  Results  from  NAVDAS-AOD 

[42]  Using  the  new  over-water  level  3  MODIS  product 
[Zhang  and  Reid,  2006]  from  January  to  May  2006,  we 
evaluated  the  performance  of  the  NRL  aerosol  data  assim¬ 
ilation  package.  In  section  6.1,  we  show  a  case  study  over 
the  west  coast  of  Africa;  in  section  6.2,  we  assess  the  impact 
of  aerosol  assimilation  on  NAAPS  performance  globally 
with  the  use  of  AERONET  and  MODIS  data. 

6.1.  An  African  Case  Study 

[43]  To  demonstrate  the  impact  of  data  assimilation  on 
NAAPS  fields,  we  use  African  dust  as  an  example.  Figure  2a 
shows  the  Terra  MODIS  true  color  image  over  the  west 
coast  of  Africa  for  30  May  2006.  A  large  dust  outbreak  is 
observed  off  the  coast  of  North  Africa.  Figure  2b  shows  the 
daily  1®  (latitude/longitude)  averaged  r  from  the  combined 
(Terra+  Aqua)  data  assimilation  quality  level  3  MODIS 
aerosol  product.  Consistent  with  Figure  2a,  heavy  dust 
plumes  with  r  around  0.5  appear  over  the  west  coast  of 
North  AfHca.  Note  that  as  mentioned  in  section  2,  stringent 
data  screening  methods  were  applied  to  the  operational  level 
2  MODIS  aerosol  product  during  the  process  of  construct¬ 
ing  a  daily  gridded  MODIS  aerosol  product  that  meets  the 
data  assimilation  standard.  Therefore,  we  expect  a  much 
reduced  spatial  coverage  from  the  newly  developed  MODIS 
level  3  aerosol  product,  compared  with  the  operational  level 
2  MODIS  aerosol  product.  Figure  2c  shows  the  r  plot  from 
the  NAAPS  nonassimilation  run  (without  use  of  NAVDAS- 
AOD)  for  1200  UTC,  30  May  2006.  The  NAAPS  non¬ 
assimilation  run  underestimated  both  the  strength  and  spa¬ 
tial  coverage  of  the  heavy  dust  plumes  shown  in  Figures  2a 
and  2b.  Figure  2d  shows  the  result  after  continuous  NAAPS 
assimilation  runs  from  1  January  2006.  As  illustrated  in 
Figure  2d,  the  NAAPS  model  run  with  the  NAVDAS-AOD 
included,  captures  both  high  r  regions  over  the  coast 
regions  and  the  transport  path  of  the  dust  plumes  across 
the  Atlantic  Ocean. 

6.2.  Global  Statistics 

[44]  To  evaluate  the  impact  of  the  data  assimilation  process 
on  NAAPS  globally,  we  compared  results  to  AERONET  data 
[Holben  et  al.,  1998].  Since  only  over-ocean  MODIS  aerosol 
retrievals  were  used  in  this  study,  our  comparisons  were 
limited  to  the  coastal  and  island  AERONET  sites.  Note  that 
AERONET  data  are  point  source  data  and  may  not  neces¬ 
sarily  represent  1  ®  x  1  ®  bin  averaged  mean  AOD  values.  A 
total  of  5  months  of  AERONET  optical  depth  data  from 
January  to  May  2006  were  included,  and  one  pair  of  NAAPS 
and  AERONET  data  were  selected  for  in  the  analysis  when 
the  time  gap  between  the  NAAPS  and  AERONET  data  was 
within  ±30  min. 

[45]  Figure  3a  shows  the  scatterplot  of  NAAPS  non¬ 
assimilation  run  and  Sun  photometer  r.  The  overall  corre¬ 
lation  between  the  two  variables  is  0.49.  The  absolute 
errors,  estimated  using  all  data  points  and  separately  using 
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Figure  2.  (a)  MODIS  true  color  image  for  30  May  2006;  (b)  daily  averaged  MODIS  (Terra  +  Aqua)  r  for 
the  same  day  as  Figure  2a.  Notice  that  with  the  stringent  data  cleaning  processes,  Figure  2b  may  be  different 
from  the  plot  that  is  created  using  the  operational  level  2  MODIS  aerosol  product;  (c)  NAAPS  r  analysis  for 
1 200  z,  30  May  2006,  without  the  use  of  NAVDAS-AOD;  (d)  NAAPS  r  analysis  for  1 200  z,  30  May  2006, 
with  5  months  of  data  assimilation. 


data  points  that  have  AERONET  r  larger  than  0.2  are  0.11 
and  0.23,  respectively.  However,  these  statistics  are  mis¬ 
leading  in  that  a  multimodal  data  distribution  pattern  is 
observed.  The  vertical  population  (Sun  photometer  r  <  0.5, 
NAAPS  r  ':$>  0.5)  represents  events  that  are  overpredicted 
by  NAAPS.  The  horizontal  population  (NAAPS  r  <  0.5, 
Sun  photometer  r  »  0.5)  represents  events  that  are  missed 
or  underpredicted  by  NAAPS.  Both  errors  may  also  be  due 
to  an  aerosol  event  in  NAAPS  that  is  phase-shifted  from  its 
real  location.  This  is  particularly  likely  near  and  downwind 
of  underreported  regions  of  the  world  such  as  Afnca  and 
South  America. 

[46]  Nevertheless,  once  data  assimilation  is  included,  we 
find  a  remarkable  improvement.  Figure  3b  shows  the 
scatterplot  of  AERONET  versus  NAAPS  r  for  the  NAAPS 
data  from  the  NAVDAS-AOD  mode  runs.  With  the  inclu¬ 
sion  of  MODIS  data,  the  correlation  between  NAAPS  and 
AERONET  r  is  improved  to  0.73,  and  more  importantly  we 
found  a  40%  reduction  in  absolute  errors  for  all  data  points 
and  for  data  points  with  AERONET  r  larger  than  0.2. 

[47]  To  increase  spatial  coverage  over  the  ocean,  we  also 
compared  NAAPS  data  against  the  level  3  MODIS  data  that 
were  developed  by  Zhang  and  Reid  [2006].  As  an  example. 
Figure  4a  shows  the  scatterplot  of  over-ocean  MODIS  r 
versus  NAAPS  r  (from  nonassimilation  run)  for  the  period 
of  March  to  May  2006  were  used  (the  NAAPS  and  MODIS 
data  from  the  full  range  of  the  study  period  is  not  shown,  as 
the  data  volume  is  enormous  and  similar  plots  like  Figure  4 


were  found  for  every  single  month).  As  with  AERONET 
comparison,  a  multimodal  pattern  is  also  observed  in 
comparing  with  daily  gridded  MODIS  data  (Figure  4a). 
This  is  due  to  the  reasons  described  in  the  previous 
paragraph,  indicating  a  need  for  improving  NAAPS  forecast 
accuracy  through  various  means,  including  the  use  of  the 
NAVDAS-AOD.  Figure  4b  shows  the  scatterplot  of  MODIS 
r  versus  NAAPS  r  (with  the  NAVDAS-AOD  included)  for 
the  same  time  period  as  Figure  4a.  Since  MODIS  data  were 
used  in  the  NAAPS  assimilation  mode  runs,  therefore, 
instead  of  comparing  MODIS  data  with  the  NAAPS  data 
from  the  assimilation  mode  runs,  we  compared  MODIS  data 
with  the  NAAPS  data  from  the  6  h  forecast  mode  runs.  The 
NAVDAS-AOD  is  used  to  modify  the  NAAPS  initial 
condition  at  t  =  0,  using  MODIS  data  from  — 3  h  <  t  < 
+3  h.  NAAPS  is  then  run  in  forecast  mode  to  t  =  +6  h,  and 
the  result  are  compared  with  MODIS  observations  having 
+3  h  <  t  <  +9  h.  With  the  use  of  the  NAVDAS-AOD,  even 
after  6  h  forecast  run,  the  correlation  is  improved  from  0.57 
to  0.81,  and  the  absolute  error  reduced  by  40-50%  for  all 
data  points  and  for  data  points  with  AERONET  r  >  0.2  as 
well.  Figure  3  and  4  indicate  that  with  careful  QA  and  QC 
processes,  satellite  data  can  benefit  aerosol  forecasts. 

[48]  Figure  4  demonstrates  the  influence  of  aerosol  data 
assimilation  through  point  comparisons  between  MODIS 
and  NAAPS  r.  Figure  5  further  illustrates  the  regional 
differences  in  NAAPS  model  runs  with  and  without  the 
inclusion  of  the  NAVDAS-AOD.  Figure  5a  shows  the 
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Figure  5b  shows  the  averaged  NAAPS  r  from  NAAPS 
nonassimilation  runs  for  the  same  period  as  Figure  5a.  The 
NAAPS  nonassimilation  runs  could  reproduce  dust  aerosol 
fronts  shown  in  Figure  5a.  However,  the  magnitudes  of  the 
dust  optical  depth  are  underestimated  for  both  regions. 
Furthermore,  the  NAAPS  nonassimilation  runs  miss  the 
major  biomass  burning  and  the  pollutant  aerosol  plums  that 
are  visible  in  Figure  5a.  Figure  5c  shows  the  three  month 
averaged  NAAPS  r  from  the  assimilation  runs.  The  NAAPS 
assimilation  runs  successfully  reproduce  both  the  spatial 
distributions  and  magnitudes  of  the  dust  episode  and  are 
also  able  to  simulate  the  smoke  and  pollutant  aerosol 
plumes  missed  by  the  NAAPS  nonassimilation  runs. 

[49]  We  also  evaluated  the  performance  of  NAAPS 
6-h  forecasts.  Figure  5d  shows  the  3 -mo nth  averaged 
NAAPS  6-h  forecast  runs  after  ingesting  observations  at 
the  analysis  time  t  =  0.  Similar  patterns  of  r  are  shown  as 
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Figure  3.  (a)  AERONET  versus  NAAPS  r  for  5-month 
(January-May  2006)  NAAPS  nonassimilation  run; 
(b)  AERONET  versus  NAAPS  r  for  5-month  (January- 
May  2006)  NAAPS  run  with  the  aerosol  data  assimilation 
process. 


3  month  averaged  (March-May  2006),  1°  x  1°  latitude/ 
longitude  gridded  MODIS  r  from  the  MODIS  level  3  prod¬ 
ucts  used  in  this  study.  Outliers  are  more  damaging  to 
NAVDAS-AOD  analysis  quality  than  data  losses.  There¬ 
fore,  in  the  newly  developed  MODIS  level  3  aerosol 
products,  stringent  data  screening  processes  (as  explained 
in  section  2)  were  applied  to  reduce  noise  and  outliers. 
Hence,  the  3  month  averaged  MODIS  r  plot  showed  in 
Figure  5a  could  be  different  from  the  plot  generated  using 
the  operational  MODIS  level  2  product.  For  example,  there 
are  data  gaps  (in  black)  over  the  high-latitude  southern 
hemisphere,  due  to  high  frequency  of  cloud  cover  over  that 
region,  which  may  not  be  observed  if  the  plot  is  generated 
using  operational  MODIS  level  2  aerosol  product.  In  Figure 
5a,  large  dust  aerosol  fronts  are  visible  off  the  West  Coast  of 
Aftica  and  East  Coast  of  Asia.  Biomass  burning  smoke  is 
observed  over  Central  America  and  the  northwest  coast  of 
North  America,  and  pollutant  aerosols  are  found  over  India. 
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Figure  4.  (a)  MODIS  versus  NAAPS  r  for  3-month 
(March  May  2006)  NAAPS  nonassimilation  run;  (b) 
MODIS  versus  NAAPS  r  for  3-month  (March  May  2006) 
NAAPS  6-h  forecast  run  (t  =  +6  h).  The  NAVDAS-AOD  is 
used  to  modify  the  NAAPS  initial  condition  at  t  =  0. 
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Figure  5.  (a)  Three-month  (March -May  2006)  averaged  MODIS  (Terra  +  Aqua)  r.  Notice  the 
newly  developed  MODIS  level  3  product  is  used.  Therefore,  there  could  be  discrepancies  between 
Figure  5a  and  the  image  that  is  created  using  the  operational  MODIS  level  2  aerosol  product  (see 
section  6.2);  (b)  3-month  (same  months  as  Figure  5a)  averaged  NAAPS  r  (total)  from  NAAPS 
nonassimilation  run;  (c)  3 -month  (same  months  as  Figure  5a)  averaged  NAAPS  r  (total)  from 
NAAPS  assimilation  run.  (d)  Three-month  (same  months  as  Figure  5a)  averaged  NAAPS  r  (total) 
from  6-h  NAAPS  forecast  run  after  assimilation  run  at  t  =  0. 


Figure  5c  with  slight  reductions  in  magnitude,  indicating 
that  the  NAVDAS-AOD  improves  model  performance  at 
even  the  6-h  forecast  timescale.  The  globally  averaged  r 
over  oceans  are  approximately  0. 1 3,  0.06,  0. 1 2,  and  0. 1 2  for 
Figures  5a- 5d,  respectively.  NAAPS  nonassimilation  runs 
underestimate  the  magnitude  of  globally  averaged  r  over 
oceans,  while  with  the  use  of  NAVDAS-AOD,  NAAPS 
aerosol  forecasting  performance  is  improved. 

[so]  As  NAAPS  is  an  operational  aerosol  forecasting 
model,  we  are  not  only  interested  in  the  improvement  in 
NAAPS’  aerosol  analysis  when  satellite  data  are  available 
but  are  also  interested  in  whether  the  use  of  NAVDAS-AOD 


at  time  t  ==  0  could  help  aerosol  forecasts  for  a  longer  period. 
Figure  6a  shows  the  absolute  error  (|Ar|,  in  absolute 
difference)  between  NAAPS  and  AERONET  r  values 
(coast  and  island  sites)  as  a  function  of  forecast  hour  in 
the  forecast  mode  run.  Data  from  March  to  May  of  2006 
were  used  in  the  analysis.  The  dashed  lines  show  the  |Ar| 
values  for  NAAPS  nonassimilation  run.  The  solid  lines 
show  the  |Ar|  values  when  the  NAAPS  run  in  the  assim¬ 
ilation  mode  before  and  at  forecast  time  0  but  not  in  any 
other  time  after  forecast  time  0.  The  black  line  with 
diamond  symbols  shows  the  |  Ar|  values  for  all  data  points 
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Figure  6.  (a)  Three-month  (March  May  2006)  averaged  differences  between  AERONET  and  NAAPS 
(nonassimilation  run)  r  as  a  function  of  number  of  hour  in  forecast  mode  run.  The  black  line  with 
diamond  symbols  is  derived  using  all  data  points.  The  black  line  with  star  symbols  is  derived  using  pairs 
of  AERONET  and  NAAPS  data  with  AERONET  r  larger  than  0.2.  Light  gray  line  shows  correlation 
between  AERONET  and  NAAPS  r  as  a  function  of  forecasting  time;  (b)  3-month  (March -May  2006) 
averaged  differences  between  MODIS  and  NAAPS  (with  data  assimilation  at  time  0)  r  as  a  function  of 
number  of  hour  in  forecast  mode  run.  The  definitions  for  solid  block,  and  light  gray  lines  are  similar  to 
that  of  Figure  6a,  except  for  MODIS  and  NAAPS  analyses. 


and  the  black  line  with  star  symbols  shows  the  |At|  values 
for  data  points  with  Sun  photometer  r  larger  than  0.2. 

[5i]  When  all  data  were  used  in  the  analysis,  the  absolute 
difference  between  AERONET  and  NAAPS  r  with  the  use  of 
NAVDAS-AOD  is  0.07  at  forecast  time  0  (A  40%  reduction 
in  I  Ar|).  As  expected,  |  Ar|  values  increase  with  time  during 
the  forecast  run.  However,  even  after  48  h  of  forecast  run,  the 
absolute  error  is  about  20%  less  than  that  without  the  use  of 
NAVDAS-AOD.  The  improvements  hold  for  larger  aerosol 
events.  For  example,  the  black  line  with  star  symbols  shows 
the  same  analysis  as  the  black  line  with  diamond  symbols  for 
data  points  with  AERONET  r  larger  than  0.2,  with  a 
reduction  in  error  of  about  20%.  The  light  gray  line  in  Figure 
6a  shows  the  correlation  between  AERONET  and  NAAPS  r 


as  a  function  of  NAAPS  forecasting  time.  After  48  h  of 
NAAPS  run  in  a  nonassimilation  mode,  the  correlation  is 
reduced  from  0.75  to  0.58. 

[52]  We  repeated  the  analysis  using  MODIS  over-ocean 
data  (Figure  6b).  The  low  |  Ar|  values  at  time  0  (solid  black 
lines)  are  due  to  MODIS  data  ingested  in  the  assimilation 
process  at  time  0.  For  up  to  48  h  in  forecast  mode  run,  the 
|At|  values  for  all  data  points  and  for  data  points  with 
AERONET  r  >  0.2  are  30  40%  lower  than  the  |  Ar|  values 
estimated  using  the  NAAPS  data  from  the  nonassimilation 
run.  Note  that  in  Figure  6b,  from  t  =  0htot  =  6h,  there  is  a 
sudden  increase  in  absolute  error.  The  sudden  increase  in 
absolute  error  could  be  attributed  to  the  fact  that  over  ocean 
satellite  data  are  ingested  at  t  =  0  h,  while  at  t  =  6  h,  large 
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Figure  7.  Spatial  correlations  of  satellite  observation  (O)  minus  NAAPS  6-h  forecast  (F)  as  a 
function  of  distance.  Averaged  (every  50  km)  spatial  correlations  of  0-F  are  shown  for  four  different 
MODIS  T  ranges:  0.0  0.1  (plus  sign),  0.1 -0.2  (star),  0. 2-0.4  (diamond),  and  0.4-0. 8  (triangle),  and 
for  all  MODIS  r  ranges  (solid  black  dots).  Vertical  bars  show  one  standard  deviation  of  the  averaged 
spatial  correlation  for  all  MODIS  r  ranges.  Black  line  shows  the  SOAR  function  (equation  (7))  fit 
with  L  set  to  200  km. 


forecasting  errors  arc  propagated  from  inland  continents  to 
global  oceans  along  the  costal  lines.  Similar  to  Figure  6a, 
the  gray  line  shows  the  correlation  between  MODIS  and 
NAAPS  T  as  a  function  of  forecasting  time.  Again,  as 
expected,  the  correlation  reduces  as  forecasting  time 
increases.  Although,  even  after  48  h  the  forecast  is 
substantially  better  with  data  assimilation.  It  worth  men¬ 
tioning  that  at  time  0,  the  correlation  between  MODIS 
and  NAAPS  r  is  near  1,  which  indicates  that  NAVDAS- 
AOD  package  properly  merges  MODIS  aerosol  optical 
depth  data  into  the  analysis  fields. 

6.3.  Innovation  Checks  and  Problems 

[53]  In  this  section,  we  evaluated  the  differences  between 
satellite  observation  (O),  and  6  h  mode  forecast  (F,  at  t  =  6  h). 
Figure  7  shows  the  spatial  correlation  ofO  minus  F  (0-F)  as  a 
function  of  distance.  The  spatial  correlation  of  0-F  for  two 
given  locations  (m,n)  is  used  to  represent  C/,(m,  n)  in  equation 
(7).  Both  satellite  and  forecasting  data  from  March  to  May 
2006  were  used.  To  calculate  spatial  correlation  of  0-F  for 
two  given  locations,  wc  require  that  there  are  at  least  30  pairs 
of  0-F  values  with  identical  time  tags,  and  we  assume  a  zero 
spatial  correlation  for  any  two  locations  separated  by  more 
than  2000  km.  Figure  7  shows  the  spatial  correlation  of  0-F 
as  a  function  of  distance  for  four  MODIS  r  ranges:  0.0  0.1 
(plus  sign),  0.1  0.2  (star),  0.2  0.4  (diamond),  and  0.4 -0.8 
(triangle),  and  for  all  MODIS  r  ranges  (solid  black  dots).  The 
vertical  bars  represent  one  standard  deviation  to  the  spatial 
correlation  values  for  all  MODIS  r  ranges.  The  black  line, 
which  is  calculated  using  equation  7  with  L  set  to  200  km,  fits 
well  with  the  solid  black  dots.  This  indicates  that  L  value  of 
200  km  is  a  reasonable  estimation  of  the  background  error 
correlation  length. 


[54]  We  also  examined  the  spatial  distributions  of  0-F, 
and  O  minus  analysis  (A,  at  t  =  0  h).  Figure  8a  shows  the  map 
of  3-month  averaged  0-A,  and  Figure  8b  shows  the  square 
root  of  variance  (STD)  of  0-A  for  March  -  May  2006.  Since  a 
high  correlation  of  above  0.98  is  constantly  observed 
between  O  and  A,  it  is  not  surprised  to  see  that  the  3  month 
mean  and  STD  of  0-A  values  are  near-zero  spatially. 

[55]  Figure  8c  and  8d  show  the  3  month  averaged  0-F,  and 
the  square  root  of  variance  of  0-F,  respectively,  for  the  same 
study  period  as  Figure  8a  and  8b.  In  Figure  8c  and  8d,  the 
global  averaged  0-F  value  is  less  than  0.01  and  the  global 
averaged  STD  value  is  on  the  order  of  0.05.  However,  large 
0-F  values  are  observed  over  west  coast  of  Africa  where 
consistent  dust  plumes  are  located.  High  0-F  values  are  also 
found  over  cast  Asia.  Given  the  fact  that  only  over-ocean 
MODIS  data  were  used  in  the  assimilation  process,  the  high 
0-F  values  arc  mostly  introduced  by  propagating  of  forecast 
errors  from  land  to  ocean.  However,  the  large  0-F  values  may 
also  be  due  to  the  following  factors: 

[56]  1.  Background  error  correlation  is  assumed  to  be 
isotropic  in  this  study.  However,  this  is  not  a  good  assump¬ 
tion  for  regions  with  sharp  aerosol  fronts. 

[57]  2.  Existence  of  forecast  biases.  A  background  bias 
correction  scheme  [e.g..  Dee  and  da  Silva,  1998;  1999],  as  a 
function  of  region  and  season,  may  need  to  be  implemented 
in  a  future  version  of  NAVDAS-AOD. 

7.  Discussion,  Conclusion,  and  Future  Work 

[58]  In  this  paper,  we  present  a  new  aerosol  data  assim¬ 
ilation  package  that  is  currently  running  in  a  pseudo- 
operational  mode  with  the  Naval  Research  Laboratory, 
Marine  Meteorology  Division’s  NAAPS  model.  A  2-D 
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Figure  8.  (a)  Global  map  of  the  three  month  averaged  satellite  observation  (O)  minus  NAAPS  analysis 
(A)  values  for  March  May  2006.  (b)  Square  root  of  variance  (STD)  of  0-A.  (c)  Similar  to  Figure  8a  but 
for  0-F.  (d)  Similar  to  Figure  8b  but  for  0-F. 


variational  technique  based  on  NAVDAS  is  used  to  inte¬ 
grate  model  output  from  NAAPS  and  satellite  observations 
from  a  newly  developed  MODIS  level  3  product.  The  new 
level  3  MODIS  aerosol  product  is  generated  from  the 
operational  MODIS  aerosol  level  2  product  through  rigor¬ 
ous  QA  processes.  An  overall  assessment  of  the  NAVDAS- 
AOD  and  directions  for  future  research  are  listed  as  follows: 

[59]  1.  With  the  use  of  satellite  observations  from 
MODIS,  especially  with  careful  data  screening  process, 
our  study  indicates  that  the  NAVDAS-AOD  not  only 
improves  the  accuracy  of  analyzed  global  aerosol  distribu¬ 
tion  at  times  when  satellite  data  are  available,  but  also 
improves  the  accuracy  of  the  estimations  in  a  forecast  mode 
for  48  h  or  more. 

[60]  2.  One  major  problem  with  this  study  is  the  three- 
dimensional  redistribution  of  aerosol  plumes  when  NAAPS 
reports  an  erroneously  low  r  value.  This  happens  when 
NAAPS  misses  an  aerosol  feature.  To  address  this  problem, 
3-D  observations  from  CALIPSO  system  are  needed  for 
future  evaluation. 

[61]  3.  As  suggested  in  section  6.3,  a  background  bias 
correction  method  [e.g..  Dee  and  da  Silva,  1998;  1999] 
needs  to  be  explored  and  possibly  used  in  a  future  version  of 
NAVDAS-AOD. 

[62]  4.  This  study  focused  on  over-ocean  aerosol  data 
assimilation.  An  over-land  aerosol  data  assimilation  process 
is  being  developed  in  a  separate  project. 

[63]  5.  The  error  covariance  matrix  for  both  background 
and  observations  are  critical  for  the  data  assimilation 
process.  In  this  study,  the  errors  in  observations  are  assumed 
to  be  uncorrelated.  The  nondiagonal  terms  of  background 
error  covariance  matrix  are  estimated  using  a  SOAR  func¬ 
tion  with  a  specified  error  correlation  length  estimated  from 
the  differences  between  MODIS  and  NAAPS  r.  However, 
as  illustrated  in  the  paper,  the  error  correlation  length  varies 
as  a  function  of  location  and  time  and  is  most  significant 


when  large  aerosol  plumes  exist.  In  future  studies,  we  plan 
to  implement  a  regional  error  correlation  length  model  and, 
eventually,  a  simplified  version  of  the  real  error  correlation 
matrix  as  computations  permit. 

[64]  6.  The  NAVDAS-AOD  is  currently  running  quasi- 
operationally  and  the  NAAPS  data  with  the  assimilation  can 
be  accessed  through  www.nrlmry.navy.mil/aerosol/  in  the 
near  future. 
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